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ABSTRACT 


Thp plate bending problems irwolving cracl s have a wide range 
o-f practical a p p 1 1 cat i on s . The finite eleiTtent method with 

adaptive ruesh refinentent such that local mesh refinement tales 
place in the rones of very liigh stress gradient gives optimum mesh 
refineiTient to the solution of such problems. A posteriori error- 
estimation and an autofnatic mesh refinement strategy pl^y vital 
roles 3n developing fully automatic adaptive analysis procedures. 
A posteriori error estimate is used in the adaptive procedure to 
assess the reliability of the finite element appro. tifTiation and to 
guide adaptive r ef i nerrient . Coupled with sup>e r con vergent patch 
recovery technique^ the adaptive analysis becomes a very powerful 
tool for the solution of a wide range of corripley plate bending 
problems. 

The other important aspect of plate bending problems is their 
finite eleiTient formulation such that thicl as well as thin plates 
can be acc ofrirrtodated in a single f orrrtula 1 1 on . The Rsissner 
Mindlin theory for the analysis of plate bending problems when 
cast with mixed forrrtulation in finite element miethod allows to 
e\plore the possibility of developing special plate bending 
elemients called 'robust elemients’ which would wort for thict 
plates as well as in the limiting case of thin plate theory. The 
major difficulty in designing such elemients was the satisfaction 
of Babusl a-Br ezz 1 stability requirements. A tnuch simpler way 
which can be applied in terms of very simple algebraic equations 
exists via the medium of patch test which is (nearly) always also 



■suf -f 1 c lenl 


A robust triangular element which satisfies all the above 
requirements has been used in the present warl. An efficient 
finite element code has been developed which requires a minimal 
inforiTiation from the user. It is seen that the element worl s 
satisfactorily and the critical Clones in plate bending cases which 
include fis<ed edges of a plate, cracl tip cone etc. are 
successively more refined. Due to the readjustment of the node 
positions, though not monotonically, the solutions tend to 
e/pected values in a few steps. 



CHAPTER I 


im'RODUCTIOM 

Plates with lateral loads find numerous applications in 
engineering. The common e 'amples being afmy tant bodie^j under 
various types of loadings, aero structures, offshore platforms, 
navy ship bodies, etc. The pi-esence of cfacls in plate structures 
under such crucial applications males the stress analysis of such 
bodies very important. But it is observed that starting from 
plates without cracls to plates with cracls, the pilate bending 
problems rarely yield to closed form solutions. This ne i_essi tates 
the use of numerical methods. For plates with cracks the ct-acl- 
tips carry very high stresses. Finite element method o+fers 
Itself as a very suitable method for all such problems. Once the 
problem is formulated for finite element analysis, it becomes 
necessary to enquire about the accuracy of the results. Also, in 
the regions of high stresses which frequently occur in the plates 
with craci-s, the mesh refinement to approidmate the euact solution 
in an optimum way becomes very important. For this the adaptive 
mesh refinement procedures are used where the magnitude of the 
error (of the desired quantity) guides the further refinenient of 
the finite element mesh. 

1.1 A BRIEF DISCUSSION ON BE3«HNO OF PLATES 

The bending properties of plate de pend greatly on its 
thicFness as compared with its other dimensions. A comprehensive 


1 



review can be found in the book by Tirroshenko et al C1959I! 

If lateral deflection of a plate is small in comparison with 
its thicJ-ness, a very satisfactory approx ifT»ate theory of bending 
of plates IS developed which is called thin plate theory. All 
stress components can be expressed by deflection <say w) of the 
plate, which is a function of the two coordinates in the plane of 
the plate. This function has to satisfy a linea'" partial 
di f -ferent xal equation, which, together with boundary conditions. 
coiTipletely defines w. Thus the solution of this equation gives all 
necessary in’f'ormat ion for calculating stresses at any point of the 
plate . 

Thiel- plate theory considers plates as a three dimensional 
probleiTi of elasticity. Here, in addition to mid plane deflections 
w(;<,y>, it IS also necessary to consider shear deformations 

and shear force components S ,S in the i<:-y plane. 

H y 

In the firesent worS- based on FEM analysis a ’robust 
triangular plate bending element is used. The elemtent is based on 
indeptendent interpolations for slopes, displacements , and shear 
forces. The name robust means that the element used in the 
f oriTiulat 1 on worl-s in the case of thicl- plates as well as in the 
litTiiting case of thin plates. 

Once the finite element results are obtained, an obvious 
question arises about their accuracy. For this reason the 


subj ect 



of error estirrates for finite elenent solutions and a consequent 


adaptive analysis in wf'icl' the appro imation is successively 
refined to achseve pr edetenruned standards of accuracy is central 
to the effective use of p.ractical, engineering, analysis. 

The present wori' uses a simple error estimator coupled with 
an automatic mesh generator (limitea to r ectangular domains). 

Two methods of refining a finite elsiiient solution er'ist. The 
first IS the simple reduction of the subdivision sice which may be 
local or global { h-r ef inenient ) . The second refinement process 
increases the order of the polirHomial trial function appr ox imiat i on 
in a predefined elerrtent subdivision ( p~ref inetrient ) . Here the first 
rr>ethod is used. The reason for this choice being its simplicity 
and also the fact that for plate bending probl ems it is very 
difficult to find out higher order elements that pass through the 
convergence requi rement s successfully. 


l.a LITEfcATUKE £URVE¥ 
1.3,1 PLATE BEMDING ® 


The 

deformation of 

a thin plate element 

1 s 

coTTipl etely 

described 

by 

its lateral 

displacement, w<x,y). 

As 

the 

governing 

differential 

equation is 

of fourth order which 

be coiTtes 

a second 


order equation in the weah f o rmulat i on , an acceptable element must 
display C continuity in the limit of mesh ref ineiTient . This means 
that w and its first derivatives must be continuous across 
intereiement boundaries. But the formulation of C* element being 



very difficult an alternative was introduced in late sixties and 
treats plates as a degenerate case of a three dimensional 
approximation with a simplified strain distribution in the 
thict^ness direction CAhamad et al . . 1968. 196^1. 

It was realised that the reduction of shear terms should 
permit the return to thin plate theory. Here, however, 
difficulties have been encountered such as ' loct- mg' ( i .e . the 
reduction of di splacements to zero) or singular mechanisms of the 
elements CPugh et al . , "*978; Hinton et al.,1979II. To elitTiinate 

loc^ ing reduced integration' concepts were introduced CPawsey et 
al.,1971; ZienViewicz et al., 19713. 

On finding that reduced integration concepts were only 
partly successful, much research worl^ followed the realisation 
that the approx imtat ion involved was identical to that used b> 
Reissner and Mindlin and that the problerri could be cast as a 
iTiixed formulation involving rotations of the normal (theta), 
shear forces (S>, and lateral displacements <w) as the primary 
variables CLee-et al., 19823. Again it was observed that no 
efficient element could be developed on the basis of mi'.ed 
formulation alone CZ lenl- lewi cr: et al . , 198B3. 

Perffaps the main difficulty stemmed frotrt the fact that 
applying the Babu5(.a“Bre2c:i stability requirements necessary for 
all mixed formulations and thus designing suitable elements 
becofCies very difficult ^;Zien^ lewicc et a).. 


198S3. 



A iTuch simpler set o-f necessary requi r e/renls was deviced for 
fTired formulation via the medium of patch test and it was shown 
CZienhiewicz et al . , 19S6I3 that it was (nearly) always also 
sufficient for ensuring stability (and hence absence of both 
Inching and singularity) The general algebraic conditions of 
Babusla and Breszi were given a simpler forrri and a conceptual 
application of th<e patch tests serves to point out the instability 
of several well l<nown f orfriuiat i ans . 

In the year 1988 a new triangular plate bending element was 
oresented. This element is based on independent interpolations 
for slopesj displacemtents and shear forces and it is shown that it 
does not suffer from an> defect cotTmion to other Mindlin plate 
elements- The element is robust and therefore suitable for use in 
commercial codes without restrictive caveats TZienliewicc et al , , 
19861 - 

a. 2. 2 PLATES WITH CRACK 

It was only after the year 1960 that the stress problem of 
crach ed plates was 8:;amined from a higher otder plate bending 
theory tah ing shear deformation into account. The earliest 
attempt CKnowles et al . , 19603 was concerned with the problerri of 
an infinite plate under uniforrri unia><ial bending far from the 
cracl and was restricted to the case of vanishingly smtall plate 
thichness. The sariie problem was studied considering the plate 
thickness to be finite CHartranft et al . , 19683. All the cases in 



these investigations were based on Reissner*s theory and led to 
the conclusion that singular stresses under bending have the same 
functional form es m the plane stress case CWilliams et al-, 
1*?573. The predictions classical theory CWilliams et al . , 
19613 are Inown to be erroneous with respect to the angular 
distribution of stress around the craci tip. As a more serious 
practical consequence there is an associated error in predicting 
the stress intensity factor whose importance in fracture mechanics 
IS well Inown. Results reported CHartranft et al . , 1968; Wang et 
al, 19683 show that the stress intensity factor varies very 
rapidly in the h/a range of 0.0 to 0.E5. In fact it was shown 
that the stress intensity factor vetsus h/a curve has infinite 
slope at h/a equal to zero CHartranft et al , ^ 19683, It is seen 
that even for h/a ratio as small as O.E the bending stress 
intensity factor is 62X greater than that given earlier CKnowles 
et al . , 19603 for h/a tending to zero. Since h'^a does not appear 
as a paraiTieter in classical theory formulation^ the inadequacy of 
this theory in handling the bending case in cracl problems becomes 
apparent. In another paper CMurthy et al , 19813 the finateness of 
the plate was tal- en into account and the analysis here was based 
on differential approach. Shear deformtation was tal-en into 
account through the use of Reissner’s theory. No restriction was 
put on the h/a ratio so that the thicl-ness effect could 
automatically be tal en into account. 



l.S 3 ADAPTIVE ANALYSIS 

Tf- e subject of error estirates for finite eleirenl solutions 
and a consequent adaptive analysis in which the approx imat ion is 
successively reamed to achieve predetermined standards of 
accuracy, is central to the effective use of finite element codes 
for practical, engineering analysis. Huch of the paoneering 
mat hemat 1 cal wori- CBabusl-a et al.,197S, -jQyoj Babusl'a, Peano et 
a]., 197S? Craig ' 19S6? Kelly et al., 19S3j Ga go et al., 

1963j needs to be translated to engineering usage. The two main 
problems encountered nave been CZi ent lewi cc: et al., 19873 

1. The cost of computations associated with error estimation and 
the difficulty of implementing such computations into an 
e'-isting code structure. 

E. The virtual impossibility of embracing a fully adaptive 
st^ucturs into an e, 'listing code structure. 

The year 1967 saw a new error estimator by Z i eni- i ewi cz-Zhu 
CZient lewi cz et al., 19873. This estimator allows the global 
energy-’ norm to be well estimated and also gives a good evaluation 
of local errors. It can be combined with a full adaptive process 
of refinement or, more simply provide guidance for mesh redesign 
which allows the user to obtain a desired accuracy. 

In a later worl- the Zi enl lewi cz-Zhu error estimator was shown 
to be effective in problerrts of plate flexure. When used in 
conjunction witl^i triongular elements and an adaptive mesh 



generator allowing a frescribed si e of elerrents to be developed 
very fast adaptive convergence for results of specified accuracy 
IS achieved CZ lenV i ewi cz et al.. 19S93. 

In another woft , presented in two parts CZieni leicz et al . , 
19931;. the issue of posteriori erior estifT‘ation was discussed. A 
general recovery technique is developed for dete rtriin mg the 
derivatives (stresses) of the finite elerrient solutions at the 
nodes. The implementation of recovery technique is simple and 
cost effective. The recovered nodal values of the derivatives with 
linear and cubic elements are supe r convergent . One order higher 
accuracy is achieved by the procedure with linear and cubic 
elements but two order higher accuracy is achieved for the 
derivatives with quadratic el ement s , 

In second part of the paper a theorem is derived showing the 
dependence of effectivity index for the Zienl iewic:c-Zhu error 
estimator on the convergence rate of the recovered solution. This 
shows that with supe r con ve rgent recovery the effectivity index 
tends asymptotically to unity. The super convergent recovery 
technique developed in the first part of the paper is used in the 
computation of the Z leni- lewi cz—Zhu error estimator to demonst rate 
accurate estimation of the exact error attainable. 

i . 3 PSiESEWT WCSSK 

Present worl is an attempt to integrate the FEM analysis of 
tyeiiding of plates with cracl- s and adaptive mesh refinement. An 



efficient eletrent (rob st triangular eletrent) has been rrade use 0 + 
to solve the thin as well as th i cl flate problerrs Also 

2 1 enl X ewx c::-Zhu erfor estimatot is used for the refinenient 
firocedure to automat i ceil re-fjne the mesh wherever the e r t- o r is 
large. 

ThuSr in the present thesis, chaptef II describes the basics 
of plate benaing and finite element formulation and chapter III 
describes the adaptive mesh refinement procedure for plate bending 
probleiTis. Chapter IV presents; results of some typical problems 
and discusses their validity. Chaptet- V summarises the 

conclusions and delineates the scope for further study. 



CHAPTER II 


PLATE BENDING FORMULATION 

a.i THE REISSNER-MINDLIW PLATE FORMULATION 

Fig, 2.1 shows a general case of plale subjected to lateral 

loading. The :< and y coordinate-> are defined to be iri the mid 

plane of the plate and the r; coordinate in normal direction witli 

corresponding component displacements u, v, w. The following 

assuniptions are made to simplify the analysis : 

(a) The strains and stf esses in the c;-~di r e ct ion are neglected 

i.e., {j; “ 0, and = 0 

tS ^ 

(b) The normals to the middle plane of the plate remain straight 
during def orrriation. 

Using the above assumptions all the displacements can be 

written in terms of mid-plane rotations & , & and lateral 

>1 y 

displacement w of the mid plane. 

Thus, 

u ~ z .<& and V - z .& (S.11 

M y 

where z is rrieasured from the middle plane and the corresponding 


strains are given as 
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be written in matri;; forffi as^ 



esses can be written as. 
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introducing linear ela?>tiL relations end integrating stresses 
in the thickness direction ps'-Trats stress resultants to be 
determineo in ternis of mad plane rotations and displacements. 
Assuming isotropic behaviour (with constant E and i-> ) in the j.-y 
plane, with the lateral shear behaviour independently defined by a 
shear modulus G. 

If the stress resultant moments are de-^ined as 
tyst 

M = z dr 

•A H 


i 

jvj = [<y z dz 
y y 

- 

t va 

M Z dz 

Ky J ay 

- tya 


(H. 


or, in the matri” form, 

M = C M !i M d 
X y xy 


these can be written as 


M = DLe 


( 2 . 4 ) 


where 
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tear resultant can be evaluated similarly as 
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the frlatrl)^ 


>c y 


more convenient to define global shear strains in terms of 
rress resultants and then 


& + ds^/^v. 

K 

0 -f Sw/0y 
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d 4 Vvr * cS 





where 


c = 


Gt 


1 


0 



(£. 8 ) 


The parameter <5 is equa] to unity in the present derisation but is 
conveniently tal- en as a value of 6/5 CZ i en (■ i ewi c 2 et al . 19S8G 
to account for non-uniform shea*" stress distribution. 


The full statefjient of the problerii involves the well I- nown 
equilibriunt resultant equations. For moirient shear equilibrium, 

-f- S = L^DL6 + S ® O < E . 9 ) 

and for lateral equilibrium, 

v’^'s f q = O (E.10) 

q being the applied load intensity norrrial to the plate rotation. 

The three sets of equations (E.9), (E.7), and (E.10) define 

the plate bending problem in its generality in terms of three sets 
of variables Sj, and w when appropriate boundary conditions are 
imposed . 

Indeed, the Kirchoff thin plate bending theory is also 
available by imposing an infinite shear modulus G or putting c=0. 
It is easy to see that in this limiting case both S and & can be 
eliminated by applying the operator V to Eq.(E.9) and 



substituting Egs 2 7 and 2 10 


T^ls giv'es 


t e well 1 n own 


bilarnonic tt’in plate equation 

( V^'l’^DLV )w + q -- 0 (£.11) 

As it 15 clear that closed -form solution of the above 
formulations is rarely possible, numerical methods are adopted. 
In Finite Element hethod it is 1 nown that for any type of element 
being used, the convergence requirements niust be satisfied. 
Convergence requirements are conditions that guarantee that e; act 
solutions will be approached as more and (Tiore elements are used to 
model the structure These requirements are listed below. 

1. The displacement field within an element rrtust be 

con t 1 nuous . 

2. The element must be able to assume a state of constant 
strain. 

In the limit of mesh refinement but not necessaril> in larger 

e 1 ement s 

3. Rigid body irsodes must be present. 

4. Elements must be compatible. 

The convergerice requirement No. 4. says that Eq , (ii.tl) 
requires C* continuity. As formulation of C* element is very 
difficult in this case, miKed formulation is adopted. 

With finite values of shear flexibility c it is easy to 
eliminate the variables S and use for solution only the & and w 



variables This indeed corresponds to the usual forirulation oi 
the frobleiT in appro intate •forrr for which such devices as reduced 
integration had to be introduced to allow approximation to 
proceed towards thin, Kirchoff limit. In what follows a full 
mi-.ed ^orm with three variables will be used, which should yield 
solution for both thicl- plates and the limiting thin pilate 

sit ua 1 1 on . 

3.2 THE MIXED FINITE ELEMENT FORMULATION AND ITS STABILITY 
CZienUewicz et al - , 1968:! 

The Finite Element formulation of the mixed ptoblem follows 
standard procedure. The problem is formulated in the three sets 
of unknowns S, and w which are apptoximated by apptopriate 
shape functions and parameters as 


& = N d 

e 

S “ N S (R, 12) 

5 

w = N w 

w 

where and N are shape functjons and S, and w are nodal 

o a V 

variables - 

Now using Galerhin method with weighting functions N^ , N^ in 

the weah form on equations (8.9) ,(2.71 and (2.10) respectively 
(using Gauss integrals on the first and last equations) the 
following equations are obtained 
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In the above 
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( 2 . 13 ) 


( 2 . 14 ) 


t = O 
s 


f =- 


r N FdT 

Jp « 

t 


J n: qdft 

C3 


and 


r n'^'cn do 

J . et Gt 


a 


is a term which becomes zero in Kirchoff thin plate theory. 

In the above T represents the two moment components imposed 
on the boundary by the traction and F the imposed shear force. 
The shape functions and N^are conveniently chosen to have 
C* continuity but can be discontinuous between elements. 



CONSTRAINED TEST RELAXED TEST INFINITE TEST 



a = 0 ft = CO a='i3/6ft = 6/ 3 


OB (9, w 2, 1 D.O.F. PER NODE (SOLID SHADE FOR RESTRAINTS) 
A A S 2 D.O.F. PER NODE (SOLID SHADE FOR RESTRAINTS) 
o • e 2 D.O.F. PER NODE (SOLID SHADE FOR RESTRAINTS) 

Figure S.E Element Patch Tests on T 6/3 a = ( n^+ n^^ ) / 'rig » 
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CONSTRAINED TEST 


RELAXED TEST 



Figure 2.3 Element Patch Tests on T 6/3 B 3 « 


INFINITE TEST 



ol = 2 ft = 3 


" ^ "w ^ ^ "s ' ^ " "s ^ 




Indt'tnl, tucft ihoite allows S to be eliniinat£?d at element level 
when P 5»f 0 and the niain problem includes only the variables & and 


w . 

The choice of the shape function and variables needs, 
however, to be such that the Babuska-Brezz i conditions are 
satisfied by the systerri and stability is ensured. The most vital 
and strictly necessary condition for non-singularity of the 
equation system (2.13) is that 
n + n 

« = 2 : 1 

n 

e 

(2.15) 


and 

ft 


n 


m 


n 

V 


2 : 1 , 


where n , n and n stand for the number of variables S, w 

respectively in any single element, element patch or indeed the 
whole problem. If either condition is violated mechanisms or 

lockiiiQ b^h^vioiir obtained • , 


3.3 THE CHOICE OF TRIANGULAR ELEMENTS: A DISCUSSION 

CZienkiewicz et al., 198S1 

It is well known that any compleK geometric domain can be 

, . ♦ .lar olpments. Also from the point of view 

discretized into triangular eiemencs. 

. thp trianoular elements can be 

of further mesh refinement, the criaujux 

•.,r,tTv/ Tn the mixed finite element 

subdivided very conveniently. m 

♦ -i r Trianole is one in which six 

formulation the simplest quadratic Triang 

4 ^ptPrmine a quadratic variation of ^ 

boundary nodes are used to determine m 

Tt^rvtp S are determined with no element 

andw, and the shear resultants S are ae 



interconnc‘clion by a linear field defined fay three internal nodes. 
This possible elertient is ternied T6/3 and in Fig. 2-2 its behaviour 
is illustrated in single element constrained and relaxed patch 
tests as well as in an infinite patch. In the relaxed test only the 
minimufTi number of degrees of freedom on the boundaries of the 
element (i.e. 3) are prescribed and condition <2.15) is satisfied. 

The infinite test is done with the assumption that two elements 
are added to a region in which all values are prescribed. The 
performance at this test deterrrdnes in general the quality of 
element behaviour, depending on the values of ot and fl achieved 
(though of course they still have to obey Eq. (2.15)). In this 
particular example a=1, indicating an incipient locking 
possibility. It is evident front the above that T6/3 is not a 
satisfactory element, but examination of the results indicates 
that satisfaction of the requirements could be achieved by 
addition of & variables in the interior of the element. This is 
accounted by creating three internal nodes or bubble functions, 
giving element T6/3 B3. The functions used here are quartics of 
the form L* L L » etc. 

4 it ft 

Fig. 2-3 repeats the three tests on this element, which in all 
cases passes successfully. 

Thus in the following finite element analysis mixed 
formulation with d ,w,S^and as variables is chosen. The 
element adopted is T6/3 B3 which has six nodes on its edges and 
three internal nodes. The eiement stiffness matrix is a matri-. of 
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Ihe order of 30 x 30 - 

Now with this element being used for the plate bending 
problems, it is an obvious question as to what extent the analysis 
is correct for the stress concentration zones. For this to be 
answered, adaptive analysis comes in the picture. The following 
chapter deals with this aspect of the problem and gives an 
implementation algorithm. 



CHAPTER III 


ADAPTIVE REFINEMENT FOR PLATE BENDING PROBLEMS 

The finite element method offers a numerically appro:<iiTiate 
solution of a mathematically posed problem. The differences 

between the exact and approximate solution, e.g., errors in 
displacenierils CZ lenkiewi cz, Taylor, 19893 

e = u - u 

‘U 

and the errors in stresses 

O' 

where Ui o refer to enact displacement and stress solutions and u, 

^ refer to approximate solutions respectively, decrease as the 
size of subdivision *h* gets smaller or as ’p’, the order of the 
polynomial in the trial function used increases. This establishes 
convergence and the acceptability of the various finite element 
forms . 


This chapter has the following objectives : 

<a) Determine the error that has occured in a particular finite 
element analysis ( a posteriori error estimate ). 

(b) Give a procedure how best to refine the approximation to 
achieve results of a given desired accuracy economically. 

In general the iteration between (a> and (b) is adaptive and 
several steps are required to achieve results that are optimal. 

The exact solution to a practical probleni is not available 
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except for sorne very siniple cases. Thus rriore accurate if not 
exact results have to be extracted from the available finite 
elenient results. 


The Zierikiewicz-Zhu error estimator CZienkiewicz et al., 
19873 has proved to be economical and effective both in evaluating 
errors simply for a given analysis and as a prelude to adaptive 
processes. The essence of the procedure is to use the difference 
between the post-processed , recovered , more accurate gradients 
(stresses! </ and those given dirctly by the finite element 
solutions O', , 


1 . e . 


■e * 

<y K 


(3.1 ) 


as a measure of the local error CZienkieicz et al . , 19923, 


The error can then be evaluated in any appropriate norm 
defined in terms of derivative errors. Most frequently the energy 
is used as it is physically more appealing, though several other 
norms can fae adopted. 

It is obvious that the success of the procedure is dependent 
on the accuracy of the recovered gradients. The global h^ 
projection and simple averaging techniques were originally 
recommended by Zienkiewicz and Zhu CZienkiewicz et al . , 19873 but 
it was soon found that such recovery techniques were inadequate 
for quadratic approximations. 

Recently a highly accurate recovery process was introduced by 
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the s<<nit» auttior s CZjf'nkieicz et al., 199Z1. The process has 
super c oriver gent propc-rties and its use in the Zienkiewi C2-Zhu 
error estimator is highly beneficial. 


3.1 EVALUATION OF ENERGY NORM FOR PLATE BENDING PROBLEMS 

The problem of thick (Rei ssner-Mindl in ) plate can be written 
in the general form as given by equations (E.9) and (S.10). With 
the constitutive relations given by equations (2.4) and (2.7) and 
defining « to be shear rigidity, the energy norm of the 
displacement or of the error can be written as 

! I ej I® = J (N - M a’o"* (M - dft + J (C - ’*'-*"* (c - dft 

a D 

VX/-V(3.2) 

\ / 

Here N and S are the finite element approximations to the 

h H 

exact M and S. As the exact values M and S are not available, 
they have to be represented by some projected values of those 
obtained from the FEM solution itself giving the a posteriori 
error estimate CZienkiewicz et al . , 19893. 

3.2 SUPERCCMMVER^NT RECOVERY PROCEDURES 

The object of the recovery of the finite element derivatives 
is to find nodal parameters^ such that the smoothened continuous 
stress fisld defined by the besis fenstions u and nodal 
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Figure 3.1 Fateh Fitting 
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paraiTieter ?. & as CZienkiewic^ et al . , 1992 D 

o- M ^ (3.3) 

is niore accurate than <5' . In the above N are the same basis 
functions as the ones used for the interpolation of the 
displacements. Here <y (with superscripts) is used as a general 
gradient term which may be stress, moment, etc. such that energy 
or its equivalent is obtained. 


In the recovery process it is assumed that the nodal values 
Of belong to a polynomial expansion of the same order p as that 
present in the basis function N and which is valid over an eleirtent 
patch (the concept of patch formation is explained in Fig. (3.1) 
surrounding the particular assembly node considered. Such a patch 
represents a union of elements containing this vertex node. This 
polynofTiial expansion is used for each corriponent of (or the 

derivatives) and 


P a 


(3.4) 


where P contains the appropriate polynomial terms and a is a set 
of unknown parameters. For one dimensional elements of order p, P 
can be written as 


P = C 1 X X* . . . x*' 3 


(3.5) 


and 


C a a ' ' a 

. 4 * • 


» • a . ^ 

P>4>4 . 


(3.6) 
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In the same way writing only the complete polynomial terms, P 
can be written for two (or three) dimensions for the appropriate 
elenient order. Thus for two dimensions and linear expansion 

P = C 1 X y 1 (3.7) 

and for quadratic 

P == C 1 X y X* xy y* 3 (3.8) 

etc . 

However , for quadrilaterals in which terms additional to the 
complete polynomial occur it appears that a slight improvement of 
results is obtained by using identical terms to those occuring in 
the shape function N. Thus for a quadrilateral 


P = C 1 K y xy 3 (3.9) 

is used and similar forms for higher expansions. 

3. 3 DISCRETE SWPERCONVERGENT RECOVERY IN AN ELEMENT PATCH 

The determination of the unknown parameters a of the 
expansion given in Eq. 3.4 is made by ensuring a least square fit 
of this to the set of superconvergent or at least high accuracy 
sampling points existing in the patch considered if such points 
are available. To do this minimization is applied on 
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A Superconvengenl 
Gauss point 
Q Nodal values 




3 node (Quadratic) elements 


Figure 


3.2 Typical Element Patches 
the Least Square Fit to 
Gauss Point Values 


in One Dimension Showing 
Sample Superconvergent 




Gauss point 

• Nodal values determined by recovery procedure 

® Patch assembly point 

Figure 3.3 Computation of Super convergent Nodal Values for 
Linear and Quadratic Quadrilateral Elements 
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F(* > 


n 


H 


y. > 


* 

( ] 
p 


.y )) 


n 

~ r ( K. , y ) - P 

%.** 4. 



(3.10) 


where coordinates of a group of sampling points , 

n = mk is the total number of sampling points and k is the 

number of sampling points on each element mi (in = 

j j 

1,S,3,...mi) of the elenrient patch CZienkiewicz at al., 19923. The 
(TiinitTiization condition of F (a ) irriplies that a satisfies 


n 

s: F' ’ 


{ x, , y. )P ( X. , y. )a 


£ P ( K, , y C K ,y. ) 

\L U K L V. 

i«»A 


(3.11) 


This can be solved in matrix form as 

a =’a"‘ b <3. IS) 

where 

fS 

A = £ p‘*‘( :< ,y )P ( x. ,y. > 

I. •^1. 1. 1- 

i»i4 

and 

n 

h - £ p'*‘(x. ,y )«y^ (X. ,y. ) 

i»*4 

The number of equations to be solved on each patch is modest and 
the recovery is performed only for each vertex node. 

The procedure is therefore quite inexpensive. It is noted 
that precisely the same matrix A occurs in the solution for each 
component of and hence only a single evaluation of this is 


necessary 
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Once the parameters a are 
values of & are calculated 
coordinates into the eKfiression 


determined the recovered nodal 


by insertion of appropriate 

. * 

for O' . 

p 


Here only the nodes inside the patch are considered. It can 

be expected that all values of o* in the dotriain of the patch are 

r> 

super convergent as it fits closely to local sampling points which 
exhibit this property and is also of the correct polynomial order. 


The procedure is explained on a one dirriensional example of 
Fig. 3.2 where linear and quadratic elements are considered. It is 
well known that superconvergence of the derivatives occurs here at 
the Gauss points shown , and the appropriate fit of linear and 
quadratic polynomials over an element patch is indicated. 

In Fig, 3.3 the sampling points for patches of quadri lateral 
elements with p = 1 and 2 are shown. Here again the derivatives at 
appropriate Gauss points are known to be super convergent . 

Obviously the procedure could be used for any order of p for 
such elemrents. In Fig. 3. 4 a similar procedure is applied to 
triangular elements of linear and quadratic types. 

For triangular shapes the existence and the locations of 
superconvergent points are still a matter which does not appear to 
have been fully explored (r.athematically despite the eat-ly 
work of Moan CMoan, 19743 suggesting the existence of optimal 





e d by r e c 


alues fo 
ts 


1 nta'>jf' a 1 1 un {>0 lilts . 


Ir. tht iin(f*r foim the average of the derivative values at 
mid sidps of adjacent elements is proved to be super conve rgent . 
Fo( fictc t 3 ca3 put poses , however , this is equivalent to stating 
ttiat the centroidal point values of the element is 
s upe r con ve r gent . 

For quadratic tringles no fully supe r conve rgent pioints appear 
to be fnown though it is suggested that some derivatives (those 
parallel to the sides) are supier convergent at two side Gauss 
points. Here , a simpler but purely e!<per imental finding that the 
central values at sides are optimal , is used. 

It is observed CZienkiewics et al., 19921! that element 
patches overlap for internal midside nodes and interior nodes in 
the element. This means that such recovered values are frequently 
evaluated from two patches. As all such values are 
5 upe r c on ve r gent , an average value is used for such nodes. The 
recovered assembly (verteK) nodal values are, however , only 
determined from a single patch. 

Thus, in this way the above formulation allows one to 
evaluate a posteriori erroi^of the finite element solution which 
is used to adaptively refine the existing mesh and thus a better 
approximation to the exact solution is obtained. 



3.4 AUTOMATIC MESHH REFINEMENT 

Automatic mesh generation and refinement is an important part 
of adaptive mesh generation codes. The initial mesh generation, 
usually coarse mesh, can be quite simple in the case of simple 
domains such as rectangles, circles etc. or complex mapping 
techniques can be used in the case of irregular regions. In the 
present analysis rectangular plate bending problems are 
considered. A simple code was developed to divide the domain into 
triangular elements and output the nodal coordinates, connectivity 
etc. which is necessary for running the finite element analysis. 

The mesh refinement part will be necessary while doing adaptive 
strategy and will be an integral part of the finite element 
analysis. When the error analysis is done and the element or 
elements to be subdivided are identified the code will 
appropr lately modify the mesh and all the necessary information. 
The particular element may be refined by dividing it into two by 
creating an extra node at the centre of one of the edges and also 
subdividing the neighbouring element also as the elements have to 
be three noded trinagles. Another way of subdividing can be 
making four smaller triangles of the original triangle to be 
subdivided by creating three nodes at the middle of the three 
edges and also subdividing the three neighbours (shown in Fig. 
3.5) into two triangles each. Several other variations of 
subdividing the elements exist and in the present analysis the 
second procedure with further illustration in Fig. 3.6 is 


followed. 
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A'..i cait !i(- i.ppri thiese procedures require several neighbourhood 
a? ray 5:. s tat It as e 1 enien t- edg es , edges— nodes, edges — elements etc. so 
as to identify and modify nodal array, element connectivity, nodal 
coordinates, boundary conditions etc. 

After the fjiesh is refined and new nodes and elerrients created, 
the nodal locations are optimised so as to get better shap>e of 
t r i ang 1 e s . 

It is needless to say that the inclusion of graphic display 
of the mesh at various stages is a necessity in the development of 
the code. 

3,5 ALGORITra>l FOR ADAPTIVE MESH REFINEMENT 

To implement the above the following algorithm is employed. 
(1) Divide the region into finite element mesh. 

<S) Obtain solution by finite element method and determine error 
norrtis in various elements. If the solution is satisfactory stop 
further refinement. 

(3) Decide the element to be subdivided. 

(4) Let k denote the triangle where the subdivision is to be done. 
This is divided into four smaller triangles, yielding new nodes at 
the mid points of the sides of the trinagle in question (Fig. 
3.6a). Triangles other than k but having new nodes on their sides 
are also divided CFig. 3.6b^d)- In this way mesh refinement 
proceeds iteratively. After each iteration mesh smoothing is done 
by adjusting the nodes such that each internal node 1 ies^^^^^^^^^ a^^^ 
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ct'iitfiide of f'olygon fornied by the surrounding nodes. 


n 



y . 


1 

n 


Tn 

1 . 

y. 

l*i *■ 


where the sunimation is done over all nodes r surrounding the i 
node. This helps in proper distribution or element sizes. 


th 


(5 > Go bad, to (2 ) . 



CHAPTER IV 


RESULTS AND DISCUSSION 

cliapler describes the results of several problems 
j 1 1 us I r c.< t i ng the effectiveness of the adaptive finite element 
method presented in the previous chapters. The plate bending 
probiem is solved starting with the triangular mesh configuration. 
Afte^ the first solution is obtained, a posteriori error estimates 
are niade on different elements of the mesh. The errors on 
different elements are compared and the element having maximum 
error is subdivided. In fact it is a ma tter of conv'enience and 
suitability as to liow many elements should be subdivided at each 
stage . 

4.1 PLATE 1HTH NO CRACK 

To validate the finite element code developed a square plate 
with all edges fixed and subject to uniformly distributed lateral 
load has been considered. The details of the f>late are as given 
in Fig. 4.1a. Due to symmetry about the center lines only a 
quarter of the plate is considered to minimize computations 
< Fig. 4. 1b). 

This ptroblem is solved with the starting iiiesh conf igut ation 
shown in Fig. 4.2a. Adaptive mesh refinement is done to obtain 
better results. Values of lateral deflection w, and moments 
M along the central line parallel to the X-axis and the Y-aKis 
are shown in Figs. 4.3 to Fig. 4.6. The analytical results are 
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al •. hftwii in the figure. 


she £* f f o f s ate larger on the boundary of the plate where both 
the {“Iges are fitted a?. is indicated in the syrrirriet r i c quarter 
cciit' 1 dated. Here in this case it is desirable to go for 
s i niu I tarieous subdivision of all the elements having energy norm 
error ratios (i.e. the ratio of energy norm error for an element 
with res.pect to the average of this error for all the elements) 
greatet than one. This is followed in the piresent case. Since 
further iterations lead to large number of degrees of freedom the 
solutions for whiclt can not be obtained by the computer system 
available, meshes only for three iterations have been shown ( Fig. 
4.2). It can be seen from the results that adaptively refined 
meshes give better approximation to the analytical values. The 
firtite element results for maximum lateral deflection are within 
37. error when compared with analytical results. Ite rat ionwise 
errors with respect to the maximium lateral deflection are also 
gi ven in Fi g . 4 .2 . 

4.2 PLATE «nTH A CENTER CRACK 

The pit'oblents related to cracks in a pilate subjected to 
loading are important. As the stresses are very high near the 
crack tip C theoretically infinite ), the element size and their 
distribution play an important role in the evaluation of stresses. 
The error determination and the consequent adaptive rrie-h 
generation play a crucial role in developing an efficient finite 
element code. In this context, a square plate with a center crack 



16 analyzed for various loading conditions. 


EDGE MOMENTS i 

Iri the case of a square plate carrying a center crack 
subjected to rtiontents applied at the edges p^arallel to the crack 
sucli that the lateral deflection is zero at these edges and the 
other two edges of the plate are free, symmetry about the two 
center lines allows the plate to be analysed for the symmetric 
quarter (Figs. 4.7a and 4.7b). The starting mesh configuration 

for this case is shown in Fig. 4.8. The moment M values are 

y 

plotted along the centre line parallel to the X-Axis. Also the 
theoretical values CSih, 19773 near the crack tip are plotted 
(Fig. 4.9). The finite element values are seen to approximate the 
theoretical crack tip moment values within an error equal to 11.7 
7. for crack tip distance of 7.93 cm. With further iterations it 
can be inferred that the solution will become more accurate. 

UNIFORMLY DISTRIBUTED LOADING i 

Here in this case the plate has one pair of parallel edges 
fi;;ed and the other two edges free. As in the earlier problem, 
this one also has symmetry about center lines and hence only a 
quarter plate is considered for finite element analysis (Figs 
4.10a and 4.10fa). The starting mesh configuration is shown in 
Fig. 4.11. Adaptive mesh refinement is carried out. It may be 
f'Ointed out that in the present analy'sis the elenient having 
ma:<imu0i error is subdivided adap>tively though several elements can 
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q = 50 N / cm 
E = 20. 6eQ6 N / 
t =0.5 cm 
^ = 0.25 
L = 24 cm 


Figure 4,1 a ^ Square Plate with all Edges Fixed and Subjected 
to Uniformly Distributed Lateral Loading 




Figure 4.ib Symmetric Quarter of the Plate for the Case 
Shown in Fig. 4.1a 







Figure 4.3 Variation of Lateral Displacement w along 

Center Line Parallel to X-Axis at Different 
Stages of Mesh Refinement for the Symmetric 
Qua r t e r in Tig. .4 i b 




Figure 4.4 Variation of Lateral Displacement w along 

Center Line Parallel to Y“Axis at Different 
Stages of Mesh Refinement for the Symmetric 


Corre5pondin»j| to Fig. 'i.Fa 
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Figure 4.5 Variation of Momont, Mx along Center Line ParalTe 
to X— Axia at Different Stages of Mesh Refinement 
for the S y mm e t r i c Q u a r t e r > n F i. g . 4 , i b 







Figure 4.7a A Square Plate with a Central Crack subjected 

to Umfrffily Distributed Bending Moment at the 
Two Edges Parallel to the Crack and the 
Remaining Two Edges Free 

M = 70 N cm / cm 
E = 20.6e06 N / cm 
t = 0.5 cm 
V = 0.25 
L = 24 cm 
a = 4.0 cm 



Figure 4.7b Symmetric Quarter of the Plate for the Case 
Shown in Fig. 4.7a 



Fi§yrt 4*^ 


Starting fitsh Configuration for the Symmetric 
Quarter Shown in- F-ig* 4. ,7 b. 


central liBRAR'i 


1 1. T., KANPUR 
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fig. 


Figurt ^ 



A Square Plate with a Central Crack Subjected 
to Unifornrily Distributed Loadingr two Edges 
Parallel to the Crack Fixed and the Remaining 
Two Edges Free 



V C/L 



K 


L / 2 



. 10 t» Symmetric Quarter of the Plate for the Case 
Shown in Fig. 4.10a 




iT 


! i i f i i 


1 ) 1 > ~r ^ ‘ ^ ^ ^ 


starting 

Quarter 


Mesh Configuration for the Symmetric 
Shown in Fig. 4.10b 


Figure 4.11 



Figurt 4.ie« FiP«t «#*h Conf 


Figure ^.12 Mesh configurations at different stages of 

Adaptive Mesh Refinement for the Symmetric 
Quarter in Fig. 4.10.b 



F-^ee 
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Figure 4.1£b Second Mesh Conflguretion 




Figurt 4.1£d Fourth Ho*h Configuration 
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Figure 4.15 Variation of Ha:;imum Energy Norm Error with Global 
Number of Degrees of Freedom for the Symmetric 
Quarter of the Plate Shown in Fig. 4.10b 
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igure 4.16 Variation of Crack Tip Moment M^ with r 
Different Stages of Mesh Refinement for 
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be subdivided at a time. This is done so as to observe the 
effectiveness of the adaptive strategy more clearly. It is seen 
that the elements near the crack tip are subdivided as the 
itf'Ccttiori pfoceed?!. Five successive mesh configurations for the 
specified boundary conditions and subjected to a. uniformlv 
distributed load are shown in Figs. 4.'12a to 4.-12e. Fig. 4.12f 
gives mesh configuration for tenth refinement. The corresponding 
plots of lateral deflection w and moment along the central line 
parallel to the X-axis are shown in Figs. 4.13 to 4.14. 

It can be seen from the figures ( Figs. 4.12a to 4.12f) that 
the mesh near the track tip is adaptively getting refined finer 
and finer. This befiaviour conforms to the expected one. The 
bending moment M when plotted along the center line parallel to 
the X-Axis can be seen to be increasing to a larger value near the 
crack tip. As the mesh is refined adaptively the crack tip value 

for li increases with some oscillatory behaviour. It is clear 

y 

from the figure (Fig. 4.14) that the results of more and more 

refined meshes will lead to better results. The variation of the 

values of moment M at the crack tip with respect to different 

y 

refinements when plotted (Fig. 4.16) shows the nature of its 
variation as refinement proceeds. It indicates that as the 
refinement is proceeding with smaller and smaller elements neat 
the crack tip the value of M at the crack tip is increasing which 
is also desirable. Also the energy norm error when plotted 
against the number of degrees of freedom in different iterations 
(Fig.4.15) it is seen that it decreases to lower values as the 



fiTH-nisMit pr o L eed-!;, IhuuQh again some oscillations are present. 



CHAPTER V 


CONCLUSIONS AND SCOPE OF FURTHER WORK 

5.1 CONCLUSIONS 

Based on the adaptive finite element plate bending analyses 
for the cases without crack and with crack the following 
conclusions can be drawn 


Based on the energy norm errors effective finite element mesh 
refinement techniques can be developed. Local super convergent 
patch recovery techniques coupled with adaptive mesh 
refinement procedures allow for local mesh refinement. 

As the mesh refinement proceeds the finite elenient results 

give better approximations to the solutions. Also 

. K* decresises with 

observed thst th© fii^Kifriurri energy norrri erro 

some oscillations as the refinement process proceeds. 


5.2 SCOPE OF FURTHER WORK 

It was observed in the course of the present 

element matrix gives very large number of *3^° 

* / 0 *p ^ 

freedom. This makes the global stiffness matri:-' 
order (of the order of 1797 x 1797 with 96 elements). 


that the 
degrees of 
very high 
Handling 


„.r, liter memory and 

such a large size mat rix require salot of co 

e beginning will 

budget. So element matrix condensation in tn 

t-ar factor was the 

allow more refinement to be carried out. Another 

element matrix. 

presence of zero diagonal elements in the 

,ij-iich can take care 

Development of a simultaneous equation solver 
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bl 


of such matrices as well as reduce the memory requirement for 
storage can also be a further investigation. 


t h e 
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